function diff = J_r(m,param)

    % calculates J in replacement region

    diff = (((1-param.omega_1).*m./rho_fun(1,'R',param) - param.omega_0./rho_fun(0,'R',param)) ...
            + param.J_rep_1*m.^(param.gamma_R_1) + param.J_rep_2*m.^(param.gamma_R_2)) ...
            +J_0_fun(m,param);

end